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ABSTRACT 

We develop a new method for estimating and removing the spectrum of the sky 
from deep spectroscopic observations; our method does not rely on simultaneous mea- 
surement of the sky spectrum with the object spectrum. The technique is based on 
the iterative subtraction of continuum estimates and Eigenvector sky models derived 
from Singular Value Decompositions (SVD) of sky spectra, and sky spectra residuals. 
Using simulated data derived from small telescope observations we demonstrate that 
the method is effective for faint objects on large telescopes. We discuss simple methods 
to combine our new technique with the simultaneous measurement of sky to obtain sky 
subtraction very near the Poisson limit. 



Subject headings: methods: data analysis — techniques: spectroscopic — galaxies: 
redshifts 
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1. Introduction 

Light from the night sky interferes with the 
analysis of the spectra of all but the bright- 
est astronomical objects. For the determination 
of galaxy redshifts this limit was reached about 
forty years ago ( [Humason et al. 1956| ). With the 
advent of digital spectrographs (e.g. |Westphal et 



al. 1975 , Shectman and Hiltner 1976|) simultane- 
ous measurement of object and night sky spectra 
became routine, thus allowing sky subtraction. 
The basic principle of sky subtraction has been 
that the sky spectrum should be observed in as 
identical a manner as possible (in time, location 
on the sky, and optical path) to the sky spec- 
trum observed in the direction of the object. As 
ever fainter objects have been observed the ac- 
curacy of the sky subtraction has become more 
critical, and increasingly sophisticated methods, 
mechanical, observational, and algorithmic, have 
been employed to match the observed sky to the 
sky in front of the object. For wide field, multi- 
object, fiber optic spectrographs this procedure 



has been especially problematic. Watson et al. 



1995 review the problem of sky subtraction with 



fiber optic spectrographs in detail. 

In this paper we present a new method for re- 
moving the influence of the night sky emission 
from the spectra of faint astronomical objects; 
this new method does not rely on the simultane- 
ous measurement of the sky with the object. 

2. Iteratively Estimating the Sky 

Removing the sky from an object + sky obser- 
vation is equivalent to estimating the sky spec- 
trum in the direction of the object at the time of 
observation. Because the sky spectrum is com- 
posed of several components the first step in the 
estimation process must be to decide which of 
these components must be accurately estimated, 
and which can be. 

Dark sky spectra may be characterized by a 
smooth continuum, plus weak absorption lines 
(from the zodiacal light) and a host of emission 



o.ao h / i 

0.15 h 

o.io i- 

0.05 j- 

o.oo IJW^v^^Xv^vjv^wA^AA^Iv 



0.30 f-( b ) 

o.io h 
o.oo \ — "-- 
-o.io t- 



4000 4500 5000 5500 6000 6500 7000 
Wavelength in angstroms 



Fig. 1. — The first two Eigenvectors in the con- 
tinuum subtracted sky model. 

lines. For many projects involving multi-fiber 
spectroscopy of faint galaxies (e.g. redshift sur- 
veys) the continuum of the object is removed be- 
fore the analysis (Kurtz and Mink 1995, hereafter 
KM98); in these instances it is not necessary to 
estimate the continuum of the sky separately; it 
suffices to estimate the continuum of the object 
+ sky. Additionally, some lines in the sky spec- 
trum are so bright that errors in their estimation, 
which are large compared with typical features in 
the object spectra, are unavoidable. 

Our new method for sky subtraction begins 
by estimating and subtracting the sky + object 
continuum; this step removes typically more than 
80% of the sky flux, as well as the object's con- 
tinuum. All through the sky subtraction process, 
we ignore pixels which correspond to the bright 
sky lines we choose not to estimate. 

Next we fit the residual sky + object spec- 
trum with a two component Eigenvector model 
(section |3]) of the continuum subtracted sky, and 
subtract this fit. Figure [l] shows the two compo- 
nent sky model used in section ^j; to first approx- 
imation the first vector represents the mean sky, 
and the second adjusts for the fluctuations in the 
OH lines. Typically we remove more than 80% 
of the flux in the residual sky spectrum with this 
step. 

The second order residual which results from 
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these steps typically shows large (wavelength) 
scale variations due to errors in the continuum 
estimation. With most of the sky lines removed 
a second order continuum fit may be made, re- 
sulting in a reflattened ( [Kurtz and Lasala 199f ) 
residual sky spectrum. 

Finally we fit a several component (in section |] 
we use ten) Eigenvector model of the reflattened 
residual sky spectrum, and subtract it. 

The residual from this last operation is the 
continuum subtracted object spectrum, along 
with any remaining noise and errors. 

3. The Eigenvector Model 



Eigenvector decomposition techniques (e.g. Courant 
and Hilbert 1924| , Press, et al. 1986Q have 
long been used to classify astronomical objects. 



Deeming 1964 used Principal Component Anal- 
ysis (PC A) to classify stars on the basis of mul- 
ticolor photometry; Brosche 1973; used PCA to 
classify galaxies on the basis of their integrated 



properties; and Kurtz 1982 used PCA to classify 
stars on the basis of their spectra. In these cases 
classification is exactly equivalent to fitting an 
Eigenvector model to the data, with the addition 
that the resulting coefficient space is partitioned 
into classes. 

More recently several authors have used PCA 
and Singular Value Decomposition (SVD, a sim- 
ilar Eigenvector technique, see the discussion in 
Press, et al. 1986 ) to classify stellar ( |Bailer-Jones 
bt al. 1998 and references therein) and galaxy 



( Promley et al. 1998 and references therein) 



spectra. 

In this work we use SVD to create two Eigen- 
vector models, first for the continuum subtracted 
sky spectra, and second for the residual reflat- 
tened sky spectra. To accomplish this task we 
begin with a set of well-observed sky spectra (we 
typically use 500 - 1000), which does not contain 
spectra with cosmic ray hits, bad focus, or other 
anomalies. We process these spectra so that each 
contains exactly the same number of counts on 



exactly the same wavelength scale. Additionally 
in all of these steps we ignore the pixels corre- 
sponding to the sky lines we choose not to esti- 
mate, (section g) Oil 5007A, NaD, Oil 6300A, 
and Oil 6363A. 

Next we subtract the continuum from each 
normalized sky spectrum using a high order (~40) 
spline fit. We decompose the continuum sub- 
tracted, normalized sky spectra into Eigenvectors 
using SVD; the first two Eigenvectors are our 
model for the continuum subtracted sky. Note 
that higher order Eigenvectors are contaminated 
by errors in the continuum fit. 

Finally we fit the two component Eigenvec- 
tor model for the continuum subtracted sky to 
our set of normalized, continuum subtracted sky 
spectra, and subtract the fit; then we reflatten 
the residual spectra using another high order 
spline fit. We then decompose the reflattened 
residual spectra using SVD; the resulting Eigen- 
vectors are our final model for the reflattened 
residual sky spectra. 

Because of the reflattening, the principal Eigen- 
vectors in this final model do not show contam- 
ination by errors in the continuum estimation; 
thus there is not an obvious cut-off point for 
the number of Eigenvectors to use in the fits. 
This cut-off will vary from project to project and 
from observing protocol to observing protocol. 
In section |3| we use ten Eigenvectors; there is a 
clear change in slope in the log-singular-value vs. 
Eigenvector-number diagram at the tenth Eigen- 
vector, indicating the transition from signal to 
noise. 

4. Testing the Method 

The purpose of our new sky subtraction method 
is to allow observations of fainter objects to 
be made with multi-object fiber-optic spectro- 
graphs, such as the 300 fiber Hectospec ([Fab- 



ricant et al. 1994, Fabricant et al. 1998a| ) on 



the converted 6.5m MMT. A typical operating 
mode for this instrument is that each half hour 
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there would be three 10 minute exposures us- 
ing the same fiber set-up; the sets of three 10 
minute spectra would be combined to remove 
cosmic rays and other shot noise. In this op- 
erating mode, with 10% of the fibers devoted to 
sky observations, several hundred defect-free sky 
spectra would be observed per night. 

This instrument is not yet in operation; here 
we therefore build sets of spectra for testing from 
the archive of spectra from the FAST spectro- 
graph ( [Fabricant et al. 1998b) ). FAST is a long 
slit spectrograph residing on the 1.5m Tillinghast 
telescope. We choose two datasets, the first of 
well observed dark skys, both to create the Eigen- 
vector models and to add to the object spectra, 
and the second a set of typical redshift survey 
spectra. 

For the sky spectra we took all sky exposures 
(exactly 1000) made between 21 Jan 1999 and 27 
Dec 1999 where the exposure time was greater 
than 500 seconds, the sun was at least 30° below 
the horizon, the moon was at least 20° below the 
horizon, and the instrumental set-up was close to 
that standardly used in redshift surveys. We re- 
moved spectra which showed signs of anomalies, 
including cosmic rays or residual noise from UV 
flooding, for a final list of 924 sky spectra. The 
mean number of counts per 1.7A pixel in these 
spectra is 569, with a factor of three variation 
about this value. 

We obtained these sky spectra from regions 
along the slit adjacent to the objects being ob- 
served; they typically measured an area on the 
sky of 3" x 30", which, for constant observing 
time, produce a flux about three times what we 
expect from a 1.5" diameter fiber on the Hec- 
tospec/MMT, after taking into account the 18.2 
times larger collecting area of the MMT. 

The galaxy spectra are the set of all spectra 
from the 15R survey ( |Geller, et al. 2000|) where 
the exposure time was greater than 500 seconds, 
the sun was at least 30° below the horizon, the 
moon was at least 20° below the horizon, and the 
redshifts were judged secure by S. Tokarz follow- 



ing manual inspection. This dataset consists of 
1648 spectra, of ~15th mag (R) galaxies, taken 
between 1994 and 1997. We then put these spec- 
tra through a fully automatic redshift reduction 
(KM98), and remove the 15 spectra where the 
fully automated result did not match the result 
confirmed by manual inspection; typically this 
difference results from unremoved cosmic ray hits 
in the spectra. We further remove 110 spectra 
where the r statistic (a measure of quality of fit, 
[Tonry and Davis 1979| ) is < 5, the redshifts from 
absorption line spectra with r > 5 are nearly al- 
ways correct (KM98); we routinely accept them 
without manual inspection. Our final test set is 
1523 typical redshift survey spectra; each original 
spectrum is reducible fully automatically. The 
mean number of counts per 1.7A pixel in these 
spectra is 277, the fluctuations about this value 
are about a factor of 1.5. 

We take each of the 1523 galaxy spectra, 
and combine it with ten different sky spectra, 
which are created by summing between 1 and 10 
randomly selected sky spectra. This procedure 
yields a sample of 15230 test spectra with sky to 
object flux ratios which vary by more than two 
orders of magnitude. 

Each of these combined spectra then has the 
sky removed by the procedure described in sec- 
tion ^, and has its redshift determined by the au- 
tomated methods of KM98. We compare these 
redshifts with those derived from the original re- 
duction of the galaxy spectra; a redshift is consid- 
ered correct if it is within 300/cm /s of the original 
value. 

Figure || shows the fraction of correct redshifts 
obtained in this test function of the flux ra- 
tio between the object and the sky, expressed in 
magnitudes. This flux ratio refers to light enter- 
ing the fiber, and is thus primarily a measure of 
the seeing convolved central surface brightness of 
the object with respect to the brightness of the 
dark sky. 

We divide the galaxy spectra into two sub- 
samples: emission line dominated spectra, and 
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Fig. 2. — The fraction of correct redshifts as a 
function of the flux ratio of the object with the 
sky, expressed in magnitudes. 

absorption line dominated spectra. The solid 
lines marked em/abs represent spectra where the 
emission/absorption line template (KM98) ob- 
tained the higher r value in the original reduc- 
tion. Emission line objects still give redshifts 
for objects about a magnitude fainter than ab- 
sorption line objects, everything else being held 
constant. The dotted line between the em/abs 
lines is the total, and represents a typical mix of 
spectra for a magnitude limited redshift survey 
at redshift 0.05 transformed to fainter apparent 
magnitude. 

The absorption line redshifts are more than 
95% correct until the objects are 2 magnitudes 
fainter than the sky and more than 50% correct 
until the objects are 3.25 magnitudes fainter than 
the sky. The emission line spectra are more than 
95% correct until the objects are 2.5 magnitudes 
fainter than the sky, and 50% correct until they 
are 4.4 magnitudes fainter than the dark sky. 

To simulate the Poisson limit we add noise to 
the original spectra corresponding to the Pois- 
son noise expected from the number of counts 
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Fig. 3. — The fraction of correct redshifts ob- 
tained for Poisson noise limited spectra, com- 
pared with the absorption line objets in figure 
H shifted 0.7 magnitudes fainter 

in the skies added to the spectra, using y/n for 
each pixel and a Gaussian random deviate gener- 
ator. This procedure simulates perfect sky sub- 
traction. We then put these spectra through the 
automated processing of KM98 and obtain red- 
shifts. 

Figure || shows the fraction with correct red- 
shifts function of brightness ratio with the 
sky. The dotted line shows the fraction cor- 
rect for the absorption line objects at the sim- 
ulated Poisson limit (i.e. with perfect sky sub- 
traction); the solid line shows the absorption line 
result from figure | shifted 0.7 mag fainter; the 
sky subtraction in this test came within 0.7 mag 
of matching the limiting magnitude imposed by 
Poisson statistics. The plots for the emission line 
objects and the sum of all objects yield very sim- 
ilar results. 

Figure |I| shows a typical absorption line spec- 
trum at several stages in the test. On top (|]a) is 
the original galaxy spectrum from the 15R sur- 
vey, its correlation with our absorption line tem- 
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Fig. 4. — A typical absorption line spectrum at 
various stages of the test; see text. 



plate gave an r value of 7.6; below it (|4]b) is that 
spectrum with the continuum subtracted. Be- 
low that (H]c) is the sum of that spectrum and a 
sky spectrum which is the sum of seven original 
sky spectra, the object is 2.2 magnitudes fainter 
than the sky. The next figure (|4]d) is the sky + 
object with the bright sky lines and the contin- 
uum removed. Second from the bottom (||]e) is 
the residual spectrum after the sky subtraction 
process, this spectrum yields the correct redshift 
with an r value of 4.2; the bottom plot (||f) shows 
the original spectrum in |||a with added random 
Poisson noise corresponding to the number of sky 
photons in ^jp, and after the continuum has been 
removed. This spectrum yields the correct red- 
shift with an r value of 6.1. 

THIS PARAGRAPH WILL NOT BE PUB- 
LISHED. Figure [5] is similar to figure but 
shows a typical emission line spectrum. The ini- 
tial correlation with the emission line template 



Fig. 5. — A typical emission line spectrum at 
various stages of the test; see text. THIS FIG- 
URE WILL NOT BE PUBLISHED. 

gave an r value of 24.6; the object is 2.5 mag- 
nitudes fainter than the sky, and the residual 
spectrum (||e) gave the correct redshift with and 
r value of 18.3, while the Poisson limited noise 
spectrum (|5]f) gave the correct redshift with an r 
value of 17.9. 

5. Discussion 

Our new method for sky subtraction can re- 
move the spectrum of the sky from faint redshift 
survey spectra with an efficiency within a fac- 
tor of two of the Poisson limit, and without any 
simultaneous measurement of the sky spectra. 

The sky spectra in the test were observed 
over one year, and contain larger variations than 
would be found in a few-night Hectospec run. 
For example the principal Eigenvector in the re- 
flattened residual model has NI 5199A as its 
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strongest feature; this line was substantially stronger 
during two months of the early summer than 
at any other time. Another Eigenvector, which 
shows a solar (zodiacal light) spectrum, has a 
clear yearly variation. Variations caused by 
changes in the instrument set-up, focus, or chip 
response (we periodically UV flood) will also af- 
fect the result. Even with these added sources 
of variance we obtain very good results for ab- 
sorption line spectra two magnitudes fainter than 
the sky, which if the Mt. Hopkins sky is V=21.46 
mag/arcsec flMassey and Foltz 200C| ) , and assum- 
ing a 1.5" diameter fiber, corresponds to an ob- 
ject with a 1.5" aperture magnitude in V of 22.8; 
to measure this would take about two hours of 
integration on the converted MMT. 

Clearly the new techniques we have described 
can be combined with simultaneous sky mea- 
surements to achieve improved results. One ap- 
proach might be to fit a continuum subtracted 
mean spectrum derived from the sky spectra in 
a single pointing (either derived from the first 
Eigenvector, or some other method of obtaining 
a robust estimate of the mean sky) , then a fit to 
a set of reflattened residual Eigenvectors, where 
the reflattened residual model is derived from a 
much larger set of reflattened residual sky spec- 
tra, perhaps those observed during an observing 
run. This procedure would, in the first pass, save 
for the continuum subtraction and the details of 
the fitting, mimic current methods for simultane- 
ous sky estimation. The second-pass fit with the 
reflattened residual model would allow subtrac- 
tion of correlated small amplitude sky flickering 
in the sky lines, and systematic effects such as 
due to bending of the fibers. 

With this combination of methods we expect 
that sky subtraction, save for the unremovable 
Poisson noise component, will not be a dominant 
limitation in deep observations with multi-object 
fiber optic spectrographs. 

We thank D. Fabricant, M. Geller, S. Tokarz, 
E. Falco, J. Roll, and W. Wyatt for discussions. 

Software implementing these techniques has 



been written in the IRAF ( [Tody 1986| ) environ- 
ment, and is available from the authors. 

REFERENCES 

Bailer- Jones, C. A. L., M. Irwin and T. von Hip- 
pel 1998. MNRAS, 298, 361 

Bromley, B. C, W. H. Press, H. Lin and R. P. 
Kirshner 1998. ApJ, 505, 25 

Brosche, P. 1973. A&A, 23, 259 

Courant, R. and Hilbert, D. 1924, Methoden 
der Mathematischen Physik I, BerlimSpringer 
Verlag 

Deeming, T. J. 1964. MNRAS, 127, 493 

Fabricant, D. G., E. H. Hertz and A. H. Szent- 
gyorgyi 1994. Proc. SPIE, 2198, 251 

Fabricant, D. G., E. N. Hertz, A. H. Szentgyor- 
gyi, R. G. Fata, J. B. Roll and J. M. Zajac 
1998. Proc. SPIE, 3355, 285 

Fabricant, D. , P. Cheimets, N. Caldwell and J. 
Geary 1998. PASP, 110, 79 

Geller, M.J., et al 2000 in preparation 

Humason, M. L., N. U. Mayall and A. R. Sandage 
1956. AJ, 61, 97 

Kurtz, M. J. 1982. Automatic spectral classifica- 
tion. Ph.D. Thesis, Dartmouth Coll., Hanover, 
NH. 

Kurtz, M. J. and D. J. Mink 1998. PASP, 110, 
934 

Kurtz, M. J. and J. Lasala 1991. in Objective 
Prism and Other Surveys, ed. A.G.D. Phillip 
and A.R. Upgren, L. Davis Press, 133 

Massey, P. and C.B. Foltz 2000, PASP, April 
2000, to appear. 

Press, W. H., Flannery, B. P. Teukolsky, S. A. 
& Vetterling, W.T. 1986, Numerical Receipes, 
Cambridge: University Press 



7 



Shectman, S. A. and W. A. Hiltner 1976. PASP, 
88, 960 



Tody, D. 1986. Proc. SPIE, 627, 733 

Tonry, J. and M. Davis 1979. AJ, 84, 1511 

Watson, F., A. R. Offer, I. J. Lewis, J. A. Bailey 
and K. Glazebrook 1998. ASP Conf. Ser. 37: 
Fiber Optics in Astronomy II, 50 

Westphal, J. A., J. Kristian and A. Sandage 
1975. ApJ, 197, L95-L98 



This 2-column preprint was prepared with the AAS M?gX 
macros v4.0. 



8 



